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LOOKING FOR O(N) NAVIER-STOKES SOLUTIONS 
ON NON-STRUCTURED MESHES 


Eric MORANO 1 
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Hampton, VA, USA 

Alain DERVIEUX 
INRIA 

BP 93, 06902 Sophia Antipolis Cedex, France 
ABSTRACT 


Multigrid methods are good candidates for the resolution of the system arising in Numerical 
Fluid Dynamics. However, the question is to know if those algorithms which are efficient for the 
Poisson equation on structured meshes will still apply well to the Euler and Navier-Stokes equations 
on unstructured meshes. The study of elliptic problems leads us to define the conditions where a Full 
Multigrid strategy has O(N) complexity. The aim of this paper is to build a comparison between the 
elliptic theory and practical CFD problems. 

First, as an introduction, we will recall some basic definitions and theorems applied to a model 
problem. The goal of this section is to point out the different properties that we need to produce 
an FMG algorithm with O(N) complexity. Then, we will show how we can apply this theory to the 
fluid dynamics equations such as Euler and Navier-Stokes equations. At last, we present some results 
which are 2nd-order accurate and some explanations about the behaviour of the FMG process. 


^his research was supported by the National Aeronautics and Space Administration under NASA Contract No. 
NAS1-19480 while the author was in residence at the Institute for Computer Applications in Science and Engineering 
(ICASE), NASA Langley Research Center, Hampton, VA 23681, by DRET Groupe 6 under contract and, INRIA and 
“Region Provence-Alpes-Cote d’Azur” (France). 
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1 INTRODUCTION 


One first important element is the mesh independent convergence speed. Hackbush, in [1] for example, 
proposes a demonstration of this property. It is done in the special case of an elliptic problem on 
structured nested meshes. We want to evaluate the properties that we must keep in order to get the 
mesh independent convergence speed when we use unstructured non-embeded meshes. 

The problem to be solved is the following: 

{ A u = f on Cl convex polygonal domain 

u Ian = 0 (i) 

u € //?( fi) and f 6 L 2 (ft) 

The discretization is a usual linear Pl-Galerkin finite element. Thus, we get a discrete space Tih 
whose dimension is equal to Nh (number of nodes), and where the subscript h indicates the mesh 
size. The resulting problem consists now in solving the linear system: 

Ah, u h = fh, ( 2 ) 

We may evaluate the discretization error thanks to the Aubin-Nitsche’s theorem and usual regularity: 

||u - < Ci h 2 ||u||fl2 < Cz h? II/Hlj (3) 

The linear system (2) may incur a lot of CPU time because of its size (large number of nodes). The 
idea is then to use a second finite element subspace Hh whose dimension Nh is less than the previous 
one (usually H = 2 h). We have then the following relationship between both spaces (also called 
grids): 

AJi 

R N " — ► R Nfl 

R T IP 

R n » — ► R Nh 

A? 

where P and R are linear interpolations (transfer operators). The iterative process can be written 
as: 

u£ +1 = Mh ul + Nh fh where: Mh = S£ (/ - P AJ? R Ah) Stf , Sh — I - uiDf' Ah 

( S defines the basic iterative smoother, v\ and v 2 the number of pre- and post-relaxations). Such 
a process converges if \\Mh\\ < 1. A very important property of this kind of method is that the 

convergence is independent of the mesh size. In order to simplify the notations (and the study) we 
rewrite Mh as the following ideal-2-grid operator [2]: 

M h = (A^ 1 - P Ajf R)(Ah St) (4) 

The norms of both factors of the right hand side of the equation (4) will determine the norm of M h : 
• The smoothing property: 

II A h Sfll < l/* 2 rj(u) , , 

lim »(«/) = 0 V ' 

is— *00 v 

depends a lot on the basic smoothing process, and, we will not give any details. 
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• The approximation property is: 


IK 1 - P A-„' R || = 0(h>) (6) 

Let us focus on (6): it takes into account the transfer operators, and overall, represents the 
difference that exists between the solution on the fine grid and the solution on the coarse grid. 

A MG scheme that exhibits these properties will result in a convergence speed that is independent of 
the mesh size: 

Vp 3i/(p) such that || Mk v h - u fc || < - u h \\ 

We may notice that demonstrating the approximation property leads to the evaluation of the following 
quantity: 

WpiPAjj'RrH-A-'W 

For nested meshes, thanks to (3), one can easily derive this from the following equality: 

PhP = PH 

On the other hand, for unnested meshes, phP is not equal to pu and this evaluation is more difficult. 
Actually, it is the same as evaluating the difference between two interpolated solutions. Zhang [3], 
thanks to Bank-Dupont’s theory, proposes such an evaluation. 

Remark : The multigrid iterative V-cycle algorithm can easily be deduced from the previous 2-grid 
algorithm recursively. It maintains the convergence speed independently of the mesh size and has an 
0(N\ogN) complexity. 

In order to apply the previous result of convergence, we propose to use a Full Multigrid (FMG) 
strategy (proposed by Brandt in [4]). A well known result is the one given by Hackbush in [1], which 
is given below: 

Theorem: We note: k the consistency order, Ci = max ( hk-i/hk) K the ratio of accuracy 

1 < k < l 

between two solutions, S the reduction factor of the MG process, i the number of MG iterations 
applied to reach the solution u\. If u * is the solution of the discrete problem, we have: 

114 - »»ll < c-c,h; with cj, = 

Assuming that C% = 2*, we deduce the exact number of cycles in each FMG phase to solve the 
lst-order problem (S' < 1/4), and the 2nd-order problem (S' < 1/8). The number of cycles i in 
each phase is constant, which leads to an algorithm that has O(N) complexity. 

Once again, the relative interpolation error [1] conditions the quality of the initialization in each 
phase. Thus, in order to stay close to the ideal scheme (where the different subspaces are nested), we 
propose to build meshes where: 

• The mesh size ratio is close to 2, 

• The triangles aspect ratio is locally comparable in the whole domain. 
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We have thus identified the different necessary ingredients to build an algorithm having O(N) com- 
plexity: 

1. A sequence of grids, . 

2. A basic smoother (ex: Jacobi, Gauss-Seidel), 

3. Intergrid transfer operators (ex: linear interpolations), 

4. MG algorithm (ex: V-cycle, W-cycle), 

5. FMG strategy. 

We may now apply it to more complex fluid dynamics problems such as the resolution of the 
N avier-Stokes / Euler equations. 

2 FLUID DYNAMICS APPLICATIONS 

We recall first the formulation of the steady Navier-Stokes equations: 

dF(W) dG{W) _ J_ fdR(W) + ?£(W)\ , W = (p,/m,pu, E) T , 
dy Re\ dx dy ) 


F(W) = 


R(W) = 



pu \ 

l 


pu 2 +p 

G{W) = 


puv 

\ (E + p)u J 




0 



T X x 


pv 
puv 
pv 


2 + p 


= ( G(W) ) 


P) 


T xy 


K dc 

are given by: , „ a x 

T “ = V ( 2 S ' If) ' Tm = I" ( 2 5y " to) ’ T -' = f, (^ + ^) 

, , U J _ „ r /*« is the Prandtl number, where p 0 , ^o, 

Re = poUoLo/no is the Reynolds number and Pv - &nd diffusivity of the considered 

and no denote respectively the characteristic densi y, y, recover the Euler equations, 

flow. It is easily seen that, if the : right-hand side is equate ^Zs^Zpose the slip condition 
The inlet conditions are defined by the farfield flow. For Euler flows, we P ^^ = and 

r / w. A A \ 


S(W) = 


0 

T x y 

T yy 


ut xv + vt. 


■jk de 
+ Trd 7 


£ 

j€K(.) 


/ 8C ./^H + /mac,-^ = 


Re 



- ] dxdy 
dy / 


( 8 ) 
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to this!’ be ‘T" two cells ’ is min ^ d by R -’ S 

farfield boundaries The 2nd order “22 e sfheT vector splitting for the 
developped by van Leer [61. We solve t r ^ “ °, bta ' ned by the use ° f ‘he MUSCL method 

[6], namely here the multistage Jacobi algorithm^, 7)2 '°" S n °°' llnear rela - x ation algorithms 

wj 0) = IVf 
For ks = 1 to nstage 

W ^" = - - cJe (b]-; w>, * .) 

*6 K(j) 

tyf?+l _ j^r( nsta 3f) 

Let us now look at the meshes. We start with an initial given fine mesh Fig.l.a. Finer meshes 


Wu = 


(9) 


are 



a. Initial 800 nodes 


b. 3114 nodes 

Figure 1: NACA0012 fine meshes. 


c. Finest 12284 nodes 


obtained by triangle subdivision (Fig.l b c) Then wp r. , . , 

[8] to build coarser meshes, from the initial one This T C ° arSemnE al *" nt,,m to Guillard 
(Fig.2.a,.b,.c). We get 6 meshes for the NACAOm^ h T “ 'f ,u “ ce of "°<le-ombeded meshes 
coarsest 19 nodes. This method all t P ro ^ e > where the finest has 12284 nodes and the 

local mesh aspectratio ^ ^ do “ t0 2 ’ a " d a “"P-able 

variables and the corrections and linear di t ° pevii ° FS t ^ are * inear interpolations, concerning the 
will be the W-cycle because it 12 ^, “ the residuals. The MO algorithm 

exist several ways to obtain 2nd-order accurate lobTons! * ' deal ' 2 " gr ' d scb<,me - Furthermore, there 

* r — on *. 

showed us that the convergence speed is hardly ind«^dCT^™tem«h h s iZ' UPW1 " d SChemeS 
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a. 223 nodes 



b. 67 nodes 



c. Coarsest 19 nodes 


Figure 2: NACA0012 coarse meshes. 


• Hemker-Koren [11] propose to get a lst-order solution with an FMG strategy and then to 
compute a certain number of DeCV-cycles: He uses the Defect Correction (DeC) algorithm 
[12], in order to solve the 2nd-order accurate following problem: Fi{W) = S. It is written: 


T x (W') = 5 , 

= f l (W N )-J r 2 {W N ) + S , N = 1,2,... 


( 10 ) 


Actually, we define the DeCV-cycling method where the lst-order problem in (10) is approxi- 
mately solved with one V-cycle. However, we do not know how many DeCV-cycles are to be 
performed and we loose the O(N) complexity. 

We propose here two different methods in order to obtain 2nd-order accurate solutions with an O(N) 
complexity algorithm [1]. 

• FMDeCV is an FMG strategy where we use DeCV-cycles in each phase, with two Jacobi sweeps 
per level (FMDeCV-2RKl). 

• FMG2 is an FMG strategy where we use on each level of the different phases the good damping 
properties of the multistage schemes (see [13]) for smoothing directly the second order accurate 
problem. 

A result of convergence of the DeCV method is given in [14] and assures that DeCV-cycling has a 
convergence speed independent of the mesh size: 

\\DeCVu a k - ul\\ < S^IK-ujH , S 2 = S\ + S\S D eC + S Dc c < 1 (U) 

where DeCVu a h is the a - th iterate of the DeCV-cycling and u 2 h is the solution of the 2nd-order 
accurate problem. 
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Remark: Desideri-Hemker in [12] show that the convergence speed SoeC of the DeC process is at 
least equal to 1/2. From (11), we should use MG lst-order accurate algorithms whose convergence 
speed Si (independent of the mesh size) is less than 1/3. Actually, the experiments showed us that 
5, >1/3 does not induce difficulty. 


3 SECOND ORDER ACCURATE RESULTS 

3.1 Euler Flows around a NACA0012 profile 





Contours interval: 0.12500E-01 

Min./ Max. Mach 0.83693,0.97262 



Contours interval: 0.50000E-01 

Min./ Max. Mach 0.54446,1.2158 



Contours interval: 0.50000E-01 

Min./ Max. Mach 0.39513,1.4854 


c. 19 nodes 


d. 67 nodes 


e. 223 nodes 


Figure 3: Euler FMG2 phases, isomach contours, M = 0.9, a = 0°. 
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Contours interval: 0.50000E-01 Contours interval: 0.50000E-01 Contours interval: 0.50000E-01 

Min./Max. Mach 0.22627,1.4630 Min./Max. Mach 0.12322,1.4793 Min./Max. Mach 0.45E-01, 1.5196 


a. 800 nodes b. 3114 nodes c. 12284 nodes 

Figure 4: Euler FMG2 phases, isomach contours, M «, = 0.9, a = 0°. 

The test-case depicted on Fig. 3 to Fig. 5 is defined by a farfield Mach number equal to 0.9, and 
a zero angle of attack. The smoother is the (4 stage) RKJ one, whose coefficients (ai = 0.14, a 2 = 
0.2939, a 3 = 0.5252, <24 = 1) are due to van Leer [15], and defines the FMG2 strategy. On Fig. 3. a, 
we present the convergence histories of the logarithm of the residual versus the number of cycles in 
each phase. The convergence is estimated as obtained when the residual decrease reaches 10 -6 . The 
first history is a 1-grid convergence on the coarsest mesh (19 nodes), the last (4-grid scheme) on the 
800 node mesh. At the end of the convergence of each phase we produce an initialization of the 
next phase by interpolating the solution on the next finer mesh. We may notice that the different 
convergence histories tend to be straight lines with the same value of slope: this allows us to say that 
the convergence speed is independent of the mesh-size and that we may use an FMG strategy. On 
Fig.3.b the residual convergence histories of the last 5 phases are depicted (the first one is a 1-grid 
convergence history with a residual decrease equal to 10 -6 ). In order to neglect oscillatory non-linear 
phenomena we choose to impose a residual decrease equal to 1 /8 (and not a defined number of cycles): 
the FMG convergence histories follow exactly the peaks of the (corresponding) phases on Fig. 3. a, and, 
the solutions are not either changed. A solution on the finest grid (Fig.l.c) is reached after only 5 
cycles (77 WU, 673 s on Convex C210 with non vectorized software). In Fig. 3 and Fig. 4 we show the 
different solutions obtained at the end of each phase: they are non-symmetrical solutions (Fig. 3), due 
to the fact that the coarse meshes are non-symmetrical. However, this phenomenon vanishes when 
the different meshes become symmetrical (Fig. 4). Another important remark is that the solution 
between the finest mesh and the next coarser one does not vary much: we may say that we get a 
nearly converged solution on the 3114 node mesh. In order to verify the previous assumption, we 
compare the FMG solution of Fig. 5 with the solution obtained after a 10 -6 residual decrease: the 
Mach number extrema are approximately the same although the isomach lines are a little bit more 
oscillatory for the FMG solution. 
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a. FMG2 


b. Converged (RES = 10 6 ) 


Figure 5: Euler FMG2 solutions, isomach contours, = 0.9, a = 0°. 


3.2 Navier-Stokes Flows around a NACA0012 profile 

The first test-case is defined by a farfield Mach number equal to 0.8, an angle of attack equal to 10 
degrees and a Reynolds number equal to 73. We use here an FMDeCV-2RKl strategy (justified by 
the mesh independent convergence of Fig.6.a). The solution (Fig.6.c) is obtained on the finest mesh 
after 4 cycles that represent 42 WU computation and a total CPU time equal to 537 s. In Fig.6.d we 
illustrate the behavior of the pressure lift, CL (drag, CD) coefficient with a solid (dash) line, versus 
the number of finest-grid iterations, up to a residual decrease on the finest grid of 10~ 12 . The points 
(cross) represent these coefficients during the FMG process, thus up to a 0.125 residual decrease. We 
can notice that the value of each of them is almost obtained at the end of the FMG phase (the error 
is equal to 64. 10 -5 for the CL coefficient and to 16.1 0~ 5 for the CD coefficient). 

The second test-case , presented in Fig.7, is defined by a farfield Mach number equal to 2, an angle 
of attack equal to 10 degrees, and a Reynolds number equal to 106. This time we use an FMG2 
strategy (more robust than FMDeCV-2RKl that needs TVD limitation and implies that the residual 
stalls from the value of 10 -3 ); this produces a solution after 7 cycles (Fig.7.c, 109 WU, 1862 s), and 
we can make the same remarks as in the previous test-cases. The next coarser mesh results in CL 
and CD values within 1% of their final values which are obtained after 7 cycles on the finest mesh 
with a related error respectively of 2. 1 0 — 6 and 2.10 -5 (Fig.7.d), 

The last test-case shows us one limitation of this method. It is defined by a farfield Mach number 
equal to 0.8, an angle of attack equal to 10 degrees and a Reynolds number equal to 500. Here again, 
we use an FMDeCV-2RKl strategy (Fig.8). We may note, once again, that the FMG convergence 
histories (Fig.8.b) look like the corresponding peaks on Fig.8. a, thus we think that the solution does 
not vary between a 10 -1 residual decrease and a 10 -6 one. However, on Fig. 9. a we note that the 
isomach lines are deformed up untill the last solution (Fig.9.c) which presents two bumps. Actually, 
since the solution differs a lot between the one obtained on the finest mesh and the other on the 
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I [Log (RES)II 



c. Density contours <*. CL and CD behaviors 


Figure 7: Navier-Stokes FMG2 solution, M «> = 2, a = 10°, Re - 106. 
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Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 

Min./Max. Mach 0. , 0.79407 Min./ Max. Mach 0. ,1.0496 Min./Max. Mach 0. , 1.0901 


c. 19 nodes d. 67 nodes e. 223 nodes 

Figure 8: Navier-Stokes FMDeCV-2RKl phases, isomach contours, M = 0.8, a = 10°, Re = 500. 
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Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 Contours interval: 0.25000E-01 

Min./ Max. Mach 0. , 1.3009 Min./Max. Mach 0. , 1.2843 Min./ Max. Mach 0. , 1.2028 

a. 800 nodes b. 3114 nodes c. 12284 nodes 


Figure 9: Navier-Stokes FMDeCV-2RKl phases, isomach contours, = 0.8, a = 10°, Re = 500. 



a. Isomach contours b. CL and CD behaviors 

Figure 10: Navier-Stokes FMDeCV-2RKl solution, Moo = 0.8, a = 10°, Re = 500. 
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next coarser mesh (Fig.9.c and .b), we may say that we did not obtain a grid-converged solution. 
Furthermore, these two humps may be justified because the meshes that we use (especially the coarse 
ones) are not adapted for such a viscous flow and may not capture the boundary layer. Our assumption 
is confirmed by the CL and CD behaviors (Fig. 10. b), where we note an error for the CL coefficient 
equal to 5.10 -2 and for the CD coefficient equal to 8.10 -3 . Moreover, we get a solution (Fig. 10. a) 
after 105 cycles (1107 WU, 14112 s), for a residual decrease of 10 -12 and which tends to confirm that 
the FMG solution is not a steady solution. 



Figure 11: Convergences 


4 CONCLUSION 

We want to point out that the use of FMG2 or FMDeCV strategy allowed us to get 2nd-order accurate 
solutions in most of cases with a limited number of operations ( 0(N ) complexity). FMDeCV-2RKl 
is more adapted to smooth problems and costs half as much as FMG2. Furthermore, we had to use an 
entropy correction technique [16] to get the above results, and occasionally a 10 -12 residual decrease 
required to increase the value of this correction to prevent the residual from stalling or to improve 
robustness on the finest grid (this increased slightly the number of cycles in each phase without 
changing the solution greatly). As depicted in Fig. 11, these types of computations do not induce any 
stall during the convergence. The main two difficulties, non-embeded meshes and the requirement of 
2nd-order accuracy, were remedied, respectively, by using a coarsening algorithm based on a Voronoi 
technique, and basic iteration techniques that were sufficient smoothers. The difficulty encountered in 
using an FMG strategy, with our meshes, increased as the Reynolds number was raised. Actually, it 
is obvious that our meshes are not adapted to these computations and that boundary layer problems 
will need the production of stretched meshes and different specialized smoothers, as suggested in [14]. 
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